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Abstract 

In the paper we propose an universal bridge functional for the closure of the Ornstein- 
Zernike (OZ) equation for the case of infinitely diluted solutions of Lennard-Jones shperes 
of different size in the Lennard-Jones fluid. Bridge functional is paprameterized using the 
data of the Molecular Dynamics (MD) simulations. We show that for all investigated systems 
the bridge functional can be efficiently papameterized with the exponential function which 
depends only on the ratio of sizes of the solute and solvent atoms. To check the parameteri- 
zation we solve the OZ equation with the closure which includes the parametrized functional 
and with the closure without the bridge functional (Hyper- netted chain closure). We show 
that introducing the bridge functional allows to obtain radial distribution functions (RDFs), 
which are close to the MD results and essentially improve predictions of the location and 
height of the first peak of the RDF. 

Keywords: Ornstein-Zernike equation, Bridge functional, Molecular Dynamics, Lennard- 
Jones fluid 

1 Introduction 

Integral Equation Theory Of Liquids (lETL) is an effective instrument for desciption of structural 
and thermodynamic properties of solutions. The main equation of lETL is the Ornstein-Zernike 
integral equation [1], which connects the direct and the total correlation functions of the particels 
of the system. Another fundamental equation which connects the direct and the total correlation 
functions is the closure relation [2]. Ornstein-Zernike equation togather with the closure relation 
allows to calculate the particle-to-particle correlation functions of the system which in turn gives an 
opportunity to calculate the main thermodynamic parameters of the system [3j . For the molecular 
systems one often uses approximations of the Ornstein-Zernike equation. The most popular of them 
are the reference interaction site model (RISM) [H [2] and three-dimensional reference interaction 
site model (3DRISM) [5l E] . There were developed efficient numerical algorithms for solving the 
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RISM and 3DRISM equations [3 El El [10]. Recently there were proposed number of methods of 
parameterizing the results of RISM and 3DR1SM calculations. These methods allow to predict 
the free energy of solvation of organic and bioactive molecules with the average error 1 kcal/mol 

im [121 [la El nails]. 

Ornstein-Zernike, RISM and 3DRISM equations give a qualitativaly correct description of the 
local stucture of the liquid. For example, RISM equations can correctly predict assymetry of 
ion solvation [T71 [18], allow to predict the stability of molecular agregates in solution [IHl 120] , 
and ligand-substrate binding [2T| [22]. However, often OZ, RISM and 3DRISM equations cannot 
reproduce many quantitative parameters such as the radial distribution functions (RDF) or the 
solvation free energy (SFE) p5l[T4l[24]. 

This is explaied by the fact, that the closure relation contains so-called bridge functional which 
is represented as an infinite sum of integrals of correlation functions and thus is practically in- 
computable [3]. Usually, to obtain the numerical solution of the equations of lETL one uses the 
empirical closures. Not all of these closures give the good coincidence with experiments. Sev- 
eral closures was obtained for the simplified models. To name a few: the Hypper Netted Chain 
(HNC) closure, where the bridge functional is simply ignored, the Percus-Yevick closure [25], the 
Martynov-Sarkisov closure [2B], the Verlet Modified closure [27] and others [28j. Parameterization 
of the molecular dynamics (MD) simulation data is another perspective method to obtaine the 
bridge functional [291 EO] . In the paper [31] it was proposed to use the closure with the repulsive 
bridge correction where the atom- atom potentials contained an additional repulsive component. 
In many cases one perform the bridge parameterization for different parts of the phase diagram, 
in particular for the density and temperature near the critical point [321 [33]. In the current paper 
we consider the bridge functional for the two-component systems with the fixed density and tem- 
perature of the system. We consider only the case of the infinitely diluted solutions, because such 
systems are the most often used for the prediction of the thermodynamic properties of biological 
molecules. In our work we consider the simplified model, where the solvent is a liquid of Lennard- 
Jones balls which have the parameter coinciding with the a parameter if the oxygen atom of 
the SPC/E water model [M] (we denote with the number 1 the sulute and with the number 2 
the solvent particles). As a solute we use the Lennard- Jones balls with diameters from 0.25(T22 to 
2(j22, where CT22 is the diameter of the solvent (here and below we call the a parameter of the LJ 
potential "the diameter"). The goal of the current paper is building the universal empirical bridge 
functional which the good coincidence with the data of the molecular dynamics (MD) simulations. 
In our work we perfrom the MD simulations for all of the investigated systems and obtain the radial 
distribution functions (RDF). By using the OZ equation we obtain the direct correlation functions 
and by using the closure relation we calculate the bridge functional for each of the systems. We 
parameterize the bridge functional with the exponential function which includes two empirical 
parameters and determine the dependency of these parameters on the size of the solvent an. In 
such a way, we obtain the universal formula for the bridge functional for the Lennard- Jones balls 
of different diameter. At the end we compare the solutions of the Ornstein-Zernike equations with 
two different closures: hypper-netted-chain (HNC) closuse which ignores the bridge functional and 
the closure which includes the proposed empirical bridge functional. We show that introducing the 
bridge functional can essentially improve the results of the calculations and make them nearlier to 
the results of the MD simulations. Out main interest is the accurate calculation of the thermody- 
namical and structural parameters of the aqueous solutions of bioactive compounds by using the 
lETL. We believe, that introducing the bridge functionals allows in many cases to substitude the 
comutationnaly expansive MD simulations with the cheaper lETL method, for example for the 
description of the interaction of osmolites with the nanoobjects and surfaces [351 ES [371 EH] and 
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Figure 1: Distibution of a Lennard- Jones potential parameters of the atoms from the 
OPLS force-filed. The data is obtained form the Gromacs 4.5.3 program packagers] (file 
oplsaa.ff/ ffnonbonded.itp) . 
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2 Desciption of the method 

2.1 Molecular dynamics simulations 
2.1.1 Description of the investigated systems 

For the molecular dynamics simulations we used the program package Gromacs 4.5.3 [13]. We 
created eight systems which contain one particle with the Lennard- Jones (LJ) potential (see the 
equation ([T])) with eight different values of an. The particle was solved with infinite dilution in 
the LJ fluid with the fixed values (J22 and €22 (these values were the same for all systems). We 
denote the solute by the digit 1 and the solvent by the digit 2. The solute and the solvent had the 
same e parameter: (en = £22)- In our work we simulated the systems with the following values of 
the (Til parameter: an = k ■ a22, k =0.25, 0.5, 0.75, 1.00, 1.25, 1.50, 1.75, 2.00. Such range of an 
parameters was used because in the popular OPLS (Optimized Potential for Liquid Simulations) 
force-field ^44j the most of a parameters lies in the range from 0.25 (T22 to 2 a22, where o"22=0.3166 
nm is the LJ parameter if the water oxygen (see Figured]). 

«,^^W^4.,((^)"-(^)=) (1) 

where i and j denote the type of particles. 

Mixed interaction parameters of the L J particles of the solute and solvent were calculated using 
the OPLS methodology: 

0"i2 = V^ii ■ 0-22 

Varaing the an we can determine the influence of the solvent at the infinite dilution in the identical 
solvent on the structural characteristics (solute-solvent radial distribution function (RDF)) and the 
shape of the bridge functional. 
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2.1.2 Thermodynamical parameters, which are used in the simulations 



In the current work we are looking for the bridge functional for the follwong solvent: LJ fluid, 
temperature T=300K, number density 33.3294 particles/nm^ (corresponds to the water density 
997.09 g/1) ,with the LJ parameters wich corresponds to the oxygen atom in the SPC/E water 
model ((722=0.316557 nm, 622=0.6500194 kJ/mol) [M]. These parameters in the reducted units 
are the following: density =1.05 73 particles/afa, T=3.8365 €22- 

2.1.3 Parameters of the molecular dynamics symulations 

We used the "leap-frog" scheme J45j for the numerical integration of the equation of motion. 
Integration step corresponded to the 0.002 ps. Cutoff radius for the LJ potential (rcutoff) was 1.2 
nm, while the LJ potential was corrected to be zero at the distance rcutoff : 



The neighbour list was created using the "cell list" method |16], which is referenced in the Gromacs 
package as the "grid search" method [15]. The neighbour list was re-newed each tenth step with 
the cutoff radius 1.4 nm. We performed the simulations in the canonical ensemble (NVT). We 
used the Berendsen termostat |17j for preserving the temperature of the system (T=300K) with 
the parameter r=2 ps. 

2.1.4 Preparing the systems and calculationg the RDF functions 

Initially we prepared the cubic cell which contain 4168 solvent LJ particles of size 5.0007 nm^ 
using the packmol program [IS]. Coordinates of the particles were optimized unsing the gradient 
descent metod [15]. After that the system was simulatated during the 1 ns in the NVT ensemble 
to establish the equilibrium. In the flnal configuration one of the solvent atoms was changed to 
the solute atom ( LJ particle with the given an). The systems with different solute atoms were 
optimized using the gradient descent method [45]. After that each of the systems was simulated 
during the 25 ns to obtain the necessary statistics. Coordinates of the particles were saved each 
0.2 ps for the further calculation of the solute-solvent RDF. RDF was calculated using the "g_rdf" 
tool from the Gromacs 4.5 package. 

2.2 Obtaining the bridge functional from the molecular dynamic sim- 
ulation data 

For the bridge functional calculation we used the RDF functions wich were obtained from the MD 
simulations. For the bridge functional calculation two types of RDF were used: solute-solvent 
RDFs {gi2{f)) and solvent-solvent RDFs ((722 ('"))• Using these functions the total correlation 
functions were calculated: hi2{r) = gi2{r) — l,z = 1,2. To calculate the bridge functional one also 
needs to know the solute-solvent direct correlation functions Ci2(r). We used the Ornstein-Zernike 
equation [H E] to obtain these functions: 



Uij{r) = 



, T <i fcutoff ) 

,r > rcutoff 



(3) 



Cl2{k) 



hi2{k) 



(4) 



l + ph22{k) 
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where p is the solvent number density (number of solvent of particles in the unit volume), Ci2(fc), 
huik), h22{k) are the Fourier-Bessel transforms of the functions Cuir), hi2{r), h22{r) correspond- 
ingly. The Fourier-Bessel transform of the total correlation functions hij{ 
formula |49l: 



is given by the following 



An 
T 



hij[r)r sm{kr)dr 



(5) 



The Fourier-Bessel transform on the discrete grid can produce artefacts which are seen as some 
oscilations of the hij{k) function near the k = 0. In the current work these artefacts were removed 
by the following method: For each pair of the neighbouring local minimum and maximum points 
/cmin and fcmax of the high-frequency oscilations the value in the point |(fcmin + ^max) was taken to 
be equal to the mean value in the extramum points: hij{^{kmin + ^max)) = |(^ij(^min) + hij{k^s.^)^ 
and afterwards the function was interpolated in the vicinity of /c = using the cubic splines basing 
on these mean values (see fig:denoise ??. 
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Figure 2: Removing the high-frequency oscilations of the functions hij{k) in the vicinity of = 0. 
In the picture is presented the function for an =£122 ■ 

The Ci2(r) can be restored from its Fourier image using the inverse Fourir-Bessel transform 



ci2(r) =^~^[ci2(A;)] 
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Ci2{k)ksin{kr)dk 



(6) 



The bridge functional -Bi2(r) can be obtained from the closure relation using the following 
formula: 

Buir) = ln^i2(r) + l3Ui2{r) - ^i2(r) + Ci2(r) + 1 (7) 

Bridge functionals for the different solute/solvent size ratios are given in FigureO Dispite the quite 
complicated shape of the bridge functionals in the first approximation they can be sufficiantly good 
approximated with the exponential function Bi2{ai,a2,r) with two fitting parameters ai and 02: 



^12(01, a2, r) = -ai exp(-a2(r - ro)) 



(8) 
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Figure 3: The bridge functionals B12 obtained from the MD simulations for the different so- 
lute/solvent size ratios (see equation ([7])) 

The distance rg is chosen is such a way that Ui2{rQ) = ISSksT, where Uu^r) is the solute-solvent 
LJ potential, ks is the Boltzmann constant, T is the temperature. The parameters Oi and 02 were 
chosen to minimize the difference between the RDF obtained from the MD simulations and the 
RDF obtained using the closure relation: 

ai,a2 = argmin||5(i2(r) - 5'ciosurc(ai, ^2, r)| | (9) 

where 5'ciosurc('^i5 0-2, = exp{— f3Ui2{r) + gi2{r) — Ci2(r) + i?i2(ai, 02, r) — 1) and the norm is defined 
by the following expression: 

00 

||5'i2(^) - 5'ciosure(ai,a2,r)|| = j (guir) - gdosmeiai, a2,r)fdr (10) 



The coefficients ai and 02 were parameterized as the functions from the solute/solvent size ratio: 
ai = ai{au/a22), 0,2 = a2(crii/o"22)- In such a way we obtain the effective bridge functional, which 
depends on one parameter - the solute/solvent size ratio. 

5l2(f^ll/f^22) = -ai(f^ll/f^22) eXp(-a2((Tii/(T22)(^ - ^o)) (H) 

To check the effectivness of the parameterization we compared the solutions of the Ornstein- 
Zernike equation with the hypper netted chain closure (i?i2(r) = 0), the solutions of the OZ 
equation with the closure which includes the empirical bridge functional ( ITTj) and the results of 
the MD simulations. The Ornstein-Zernike equations were solved using the iterative algorithm. 
To do this the functions 712 (^) = ^12 (^) ~ Ci2(r) were introduced. The following algorithm was 
used: the first approximation was chosen to be 712^^") = 0. The (n + l)-st approximation '^\2'^\r) 
was obtained from the n-th approximation using the following algorithm[7]: 

• Step 1: ci2(r) = exp(-/3f/i2(r) + 7i^^(r) + 5i2(r)) - 7i2(r) - 1; 

• Step 2: Ci2{k) = J'[ci2(r)] 

• Step 3: ^[2^^\k) = pcuik) ■ h22{k) 
. Step 4: 7irHr) = -F-M7ir^(A:)] 



6 



3 
2 



o 


-1 

0.5 1 , 1.5 2 

Figure 4: Dependence of the In(ai) on the solute/solvent size ratio 
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Figure 5: Dependence of the 02 on the solute/solvent size ratio 

3 Results 

We performed the MD calculations for the infinitely diluted solutions of the LJ balls of dimeter 
(Til in the LJ fluid which consists of the balls of diameter 022- The following solute/solvent size 
ratios were used: crii/(T22 = 0.25, 0.5, 0.75, 1, 1.25, 1.5, 1.75, 2. From the MD simulation results 
the bridge functional presented on the Figure [3] were obtained by using the expressions (jl])- ([71) 
. For each ratio crii/cr22 the obtained bridge functional was fitted with the function ([8]) where 
coefficients 01,02 were obtained as a result of solving the minimization problem ([9]). 

In Figure H] and Figure [5] one can see the dependencies of the coefficients a\ and 02 on the 
c"ii/c"22 ratio. As one can see in Figure coefficient 02 weakly depend on the 0"ii/cr22 (standard 
deviation from the mean value is 3.3 nm~^, while the mean value is 30 nm~^ ). The coefficient 
02 can be approximated by its mean value: 02 ~ 30 nm Coefficient a\ have an exponential 
dependence from the parameter 0"ii/cr22 and thus In(ai) linearly depend on o\\jo22- In Figure H] 
for all crii/(j22 ratios from 0.5 to 2 the clear linear dependence of In(ai) on o"ii/(T22 is seen: 

In(ai) = C • aii/a22 (12) 




a2 = 30±3.3nm ^ 
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Figure 6: Comparison of the RDFs obtained from the Ornstein-Zernike equation with the closure 
without bridge functional and with the bridge empirical functional with the MD simulation results. 

'7ll/o"22 = 1-75 

The exception is the ratio (Tii/(T22 = 0.25, which we do not use for the parameterization. For 
the values pf Oi, which correspond to the cTii/cr22 = 0.5, 0.75, 1, 1.25, 1.5, 1.75, 2 using the 
latest squares method the value of the C coefficient in the equation eq:log(al)=Csigma (??)as 
determined: C = 1.1754. In such a way, the final formula of the empirical bridge functional is the 
following: 

512(0-11/0-22) = -exp{-ai{r - vq) + Ccrii/0-22) (13) 

where ai=30 nm-\ C = 1.1754, Uuiro) = U.SksT. 

To check the empirical functional (|T3l) for each of the ratios 0-11/0-22 we solve the Ornstein- 
Zernike equation with two different closures: hyper-netted chain, where the bridge functional 
Bi2{r) = 0, and with the closure, which includes the empirical functional built using the formula 
([T3l) (the h22{r) was taken from the MD simulation of the pure solvent). In Figure [6] the RDF 
obtained as a result of solving the OZ equation with the closure which includes the empirical 
bridge functional ffT5]) . the RDF obtained as a result of solving the OZ equation with the closure 
without the bridge functional and the results of the MD simulation are shown. The solute/solvent 
size ratio is 0-11/0-22 = 1.75. One can see, that the RDF, obtained without the bridge functional 
predicts incorrectly the position of the first peak of the RDF and also overestimate the height 
of the peak, while the RDF, obtained with the empirical bridge functional predicts correctly the 
position of the peak and more accurately predicts its height. In Tabled] the comparison of the first 
peak height obtained by using the different closures (with the bridge functional and without it) 
with the MD results is presented. One can see that both of the closures overestimate the height 
of the first peak of RDF, but in averae the error of the closure with the bridge functional is 0.109, 
which is approximately 2.5 times lower than the error of the closure without the bridge. In Table 
[2] the comparison of the first peak position of RDF, obtained with the two different closures (with 
the bridge functional and without it) with the MD results is presented. One can see, that the 
closure without the bridge functional gives the systematic error of the first peak of RDF position 
around 0.01 nm (O.lA), while the closure, which includes the bridge functional, determines the first 
peak position with the error which is 10 times lower (0.001 nm or O.OIA). This error is negligable, 
because it is 10 times lower than the grid size which was used for the discretization of RDFs 
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Table 1: Comparison of the heights of the first peak of RDFs which were obtained by using the 
different closures with the MD simulations. 
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Table 2: Comparison of the positions of the first peak of RDFs which was obtained using the 
different closures with the MD simulations. 
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obtained from the MD simulations. Additionally one need to stress the importance of the accurate 
approximations of the first peak position of the RDF for the sovation free energy calculations, 
because for that application small changes in the first peak position lead to the essential errors in 
the calculations. 

4 Conclusions 

In the paper the new universal bridge functional for the infinitely diluted solutions of the LJ spheres 
of different size in the LJ fluid is proposed. The investigated systems were the simplified model 
of the important case of aqueous solutions of different bioactive compounds, cr-parameter of the 
solvent LJ particles was chosen to be equal to cr parameter of the water oxygen in the SPC/E water 
model. The ratios of the solute/solvent particle sizes o"ii/(T22 was variaing in the range from 0.25 
to 2 which approximetely corresponds to the distribution of the a parameters of different atoms 
in the OPLS force-field. 

For the mentioned above systems the MD simulations were performed. Using the RDFs, ob- 
tained from the MD simulations, and usign the Ornstein-Zernike equation, the bridge functionals 
were obtained, wich in turn were fitted by the exponential function (|H]), which depends on two 
parameters a\ and 02. 

Dependency of the parameters and 02 on the crii/cr22 ratio was investigated. It was shown 
that the 02 parameter weakly depends on the size of the solute, and for all of the investigated 
systems it can be put to be equal = 30 nm~^. It was shown that In(ai) correlates to the 
parameter o"ii/(T22, and it can be estimated with a good accuracy using the formula ( |T2|) . The 
final empirical bridge functional ( ITSj) was determined. This functional depends only on the 
solute/solvent particle size ratio. 

Using the iterative algorithm the solutions of the Ornstein-Zernike equation for the closure with 
the empirical bridge functional were obtained. It was shown that introducing the bridge functional 
into the closue relation reduces by 2.5 times the error in the calculation of the first peak height 
of the RDF and, in contrast to the closure without the bridge, predicts accurately the position of 
the first peak of the RDF (with the error which is 10 times smaller than the grid discretization 
step). We have shown that introducing the empirical bridge functional into the closure relation 
can essentially improve the accuracy the predictions of the first peak of RDF position. 
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